Complexity of Polarized Spatial Patterns in Large Area Square VCSEL 
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We consider pattern selection process in a wide aperture VCSEL near threshold. We show that 
for a square geometry of the laser aperture, the patterns formed at lasing threshold can be very 
complicated because of a possible misalignment between directions of an intrinsic spatial anisotropy 
of VCSEL and lateral boundaries of its aperture. The analogy with quantum billiard structures is 
established, and fingerprints of wave chaos are found. Influence of localized inhomogeneous in the 
pump current is also considered. 

PACS numbers: 42.60.Jf, 42.65.Sf 



I. INTRODUCTION 

In the last decades, semiconductor laser devices have 
played an increasing role in scientific research and appli- 
cations. Among them, vertical cavity surface emitting 
lasers (VCSELs) can be distinguished, since these lasers 
emit normal to the wafer surface in the direction of epi- 
taxial growth. One of the features of VCSELs design is 
a possibility to mount extremely large (up to hundreds 
jim in diameter) aperture with high level of spatial ho- 
mogeneity. Spatial mode structures in such wide- aper- 
ture VCSELs have been a subject of many researches 

flaiiiiaiiEini. 

In general, the mechanism of pattern formation in VC- 
SEL is very complicated and involves both a complex 
structure of a VCSEL cavity (including Bragg reflectors) 
[Tol . [Til . IT3 |. as well as peculiarities of light-matter in- 
teraction in active quantum well semiconductor layers 
[H, m, [13, m, [H However, near threshold one 

can invoke the perturbation theory and obtain normal 
forms governing the evolution of the system. It was 
shown recently, that due to a slight spatial anisotropy 
of a VCSEL cavity only a few spatial modes come into 
play at threshold and the shape of these modes can be 
analyzed by the linear approximation of the normal form 
[Tol . [Til . [T^ l (in contrast to spatially isotropic systems, 
where the whole degenerate family of modes have the 
same critical growth rate at threshold, and the selection 
process requires consideration of nonlinear competition 
even at threshold [U H^j). The investigations of the 
problem using this point of view was started in [l2| for 
VCSELs with circular aperture. In that work, transition 
from the 'flower-like' modes dictated by circular bound- 
aries to 'stripe-like' ones, which are required by spatial 



anisotropy of a VCSEL, was described. 

In [T^l, the description was restricted to the light lin- 
ear polarized in direction coinciding with one of axes 
of the intrinsic laser anisotropy. In the present article, 
we extend this approach to the more general "vectorial" 
case when the laser anisotropy is not so strong and both 
orthogonally polarized components of the field must be 
taken into account. This extension allows us to investi- 
gate the interaction of two polarization degrees of free- 
dom, that is important for the consideration of competi- 
tion of different mechanisms affecting pattern formation. 

The current research is motivated by the recent ex- 
perimental investigations [T3, [3, [3| , where it is shown 
that many aspects of spatial structures in VCSEL such 
as their transverse shape or frequency-versus-lengthscale 
dependence can be well described in a linear approxi- 
mation, which validates the 'linearized normal form' ap- 
proach, used in this article. 

We show here that this makes the structures at the 
laser threshold very complicated even in the case of a 
simple square aperture and perfectly manufactured VC- 
SEL cavity. The similar complexity was observed in the 
resent experimental investigations [1, IH, H^] , where the 
scarred structures in the square broad-area VCSEL were 
obtained and treated as coherent states of quantum bil- 
liard. So, it is possible to speak in this context about 
'quantum chaos' in a wide aperture VCSEL. 

The 'quantum chaos' (or 'wave chaos') approach, in- 
cluding the investigation of level (or eigenfrequency) dis- 
tribution in a domain with 'integrable' or 'nonintegrable' 
boundary conditions, was applied for investigations of 
microwave billiar ds l,2 5j, 26,, 27| , as well as semiconductor 
microdisk lasers [28l |29| . Very recently, the wave chaos 
was observed in a VCSEL with deformed circular cav- 
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ity geometry [30|. Though the laser emission far from 
threshold is strongly affected by nonlinearities, the opti- 
cal spectra obey the same relations as a quantum billiard, 
governed by a linear Schrodinger equation. 

We develop this approach further by considering the 
wide-aperture VCSEL near threshold as a kind of quan- 
tum billiard where the role of Hamilton operator is played 
by the operator of linearized order parameter equation 
for VCSEL. We show, that this operator is considerably 
more complicated than the transverse Laplace operator 
which is often used for consideration of quantum billiard 
systems. This leads to appearance of quantum chaos-like 
features, such as scarred orbits and nongaussian statistics 
of the levels even in the simple square geometry without 
a deformation of boundary conditions, that usually used 
for demonstrating quantum chaos features. 

The structure of the article is following: in the sec- 
tion [IT] we extend the approach used in to take into 
account both orthogonally polarized components of the 
light field. In the section UIIl we consider the pattern for- 
mation in VCSEL with a square geometry and discuss the 
influence of the polarization axes rotation against lateral 
boundaries of the cavity, as well as of slight spatial per- 
turbations of the pumping current. In the section HVl we 
interpret our results in the 'quantum chaos' framework. 



II. THE VECTORIAL EIGENVALUE 
PROBLEM 

As has been pointed in the Introduction, the eigen- 
modes and their decay (or growth) rates can be directly 
obtained from a linear stability analysis of the nonlas- 
ing zero solution of the nonlinear equations governing 
VCSEL's dynamics [l^, [O, [13] ■ I^i comparison with a 
method when the field is decomposed into transverse 
eigenfunctions of the empty cavity |3l|, [s^, this ap- 
proach takes into account properties of a gain medium 
and of Bragg mirrors composing the cavity by the strick 
way. 

For the infinitely large transverse-area devices the re- 
sulting eigenmodes are plane tilted waves (transverse 
Fourier modes) with the decay rates depending on the de- 
tuning of the longitudinal cavity resonance from the gain 
maximum and on a polarization of modes (lol . [ill . [33| . In 
contrast, the corresponding modes for the finite devices 
may be very complicated even in the scalar case. For 
example, for the circular aperture the modes can range 
from Bessel-like modes to the modes resembling more 
tilted waves [T^ . 

The basic nonlinear equations for the vectorial case 
were obtained in [ll| . The short description of these 
equations is given in Appendix]^ The linearization pro- 
cedure is a generalization of one introduced in [12| to 
vectorial case, and is described in Appendix IB] As in the 
scalar case, we linearize a complex nonlocal nonlinear 
operator, obtaining the linear but nonlocal pseudodiffer- 
ential operator (or speaking more precisely, certain eigen- 



value problem for such operator) governing the field at 
threshold. Then, we approximate the operator describing 
the action of Bragg reflectors in transverse Fourier space 
and obtain the operator containing only partial deriva- 
tives up to fourth order. The total operator obtained af- 
ter such procedure includes main peculiarities of pattern 
formation such as the dependence on detuning between 
longitudinal resonance of the cavity and gain peak fre- 
quency, and on the anisotropy of Bragg reflectors. This 
operator has the following form in the (x, y) space: 



= ^ a.ij-—-—+goa5fi{x,y) + il5n{x,y). (1) 

i,j=o ^ y 

In contrast to the scalar case, for the vectorial case 
operator O is a matrix, acting on the vectorial field 
E = (ExjEy). Therefore, the coefficients also 
2x2 matrices. We also consider inhomogeneities in the 
pumping current Sfi and in the index Sn, which play the 
role of disturbances. In the following, we mainly restrict 
ourself to only current inhomogeneities, because they are 
easier to be controlled. 

Eigenfunctions Eg and eigenvalues Xg of this operator 
are obtained by solving an eigenvalue problem: 

d-Eg{x,y)-Xg-Eg{x,y) = 0. (2) 

The eigenfunction corresponding to the eigenvalue with 
the largest real part determines a spatial field distribu- 
tion with maximal growth rate at threshold appearing 
after onset of generation. Aperture of the device has 
been simulated by zero field conditions on the bound- 
aries. Besides, because the operator O is of fourth order, 
one should introduce an additional boundary condition 
which contains spatial derivatives of the field (see Ap- 
pendix [BJ to obtain the well posed eigenvalue problem. 
This second boundary condition can be chosen by differ- 
ent wa ys, but all them leads to approximately the same 
result nj. Therefore, we select the second boundary 
condition appropriately, to simplify the solution (see Ap- 
pendix [C| . 

The eigenvalue problem Q is a generalization of more 
conventional one with transverse Laplace operator ([6]). 
The important difference however, is that the operator 
P]) is anisotropic, with a principal directions defined by 
polarization anisotropy of VCSEL cavity. It should be 
noted that x-component of the field and y-component 
of the field Ey enters to Eq. ^ as just two components 
of a single vectorial eigenfunction E^. In the following, 
we suppose that the x component of the field is always 
strongest one due to the intrinsic polarization anisotropy. 
However, the orthogonal component can not vanish be- 
cause of the mixing of both components during propaga- 
tion in the cavity. 
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III. COMPLEX SPATIAL STRUCTURES IN 
SQUARE-SHAPED VCSEL 

A. Influence of alignment of boundaries and 
ansotropy direction 

To elucidate the mechanism governing the pattern se- 
lection, we first consider the case of homogeneous pump- 
ing and index profile {5fi = 0, Sn = 0). The presence 
of the aperture is modeled by zero boundary conditions. 
In this subsection we show, that the resulting spatio- 
temporal distribution strongly depends on the alignment 
of the boundary conditions and the main anisitropy axes 
of the device. 

When boundaries and anisotropy directions are com- 
pletely aligned, the a;-component of the first mode at 
threshold is a standing wave with stripes parallel to the 
direction of its polarization (see Fig. [T]Ja), (b)). Though 
the y-component is not zero due to the polarization mix- 
ing effect by Bragg reflectors, the coupling between x and 
y polarizations is very week since for this case the Fourier 
transform of homogeneous part of O, 

4 

Oik) = E {-iy^'a^lKH■ (3) 

which is {kx, fcy)-dependent 2x2 matrix, have nearly zero 
diagonal elements on the line ky = Q and kx = 0. Thus, 
y-component is weaker in 100 times approximately for 
taken parameters. 

If the boundaries are rotated with respect to the po- 
larization anisotropy direction, the stripes remain for a 
small rotation angle a, and their direction coincides with 
the boundaries rather then with anisotropy direction, i.e. 
they are rotated with the boundaries (see Fig. [He)-(h)). 
Structures of both polarized components are similar and 
their intensity become comparable. 

However, as a reaches some critical value, the corre- 
sponding structure becomes more complicated. Thus, for 
the parameter of Fig. [T](i)-(1), the spatial distribution of 
the x-polarized component of the field resembles stripes 
in the middle of the aperture, whereas near boundaries 
the pattern is more 'square-like' one. At the same time, 
the structure of y-polarized component is more compli- 
cated and less regular. The critical angle for which stripes 
are giving place to complicated structures decreases with 
the size of the device, as well as with a value of anisotropy 
7a. However, stripes take place for relatively large size 
(composing of several tens of oscillation of the intensity) . 

It should be noted that the structures presented in 
Fig. [IJi)-(l) arc the combination of many transverse 
Fourier harmonics of different directions (sec Fig.[IJj),(l)) 
in contrast to two travelling waves composing stripes for 
small a (Fig. [Hb),(f)).We show in the next section by 
considering the corresponding eigenvalue statistics, that 
such structures demonstrate some features which make 
them close to quantum chaos wavefunctions . 
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FIG. 1: The structures at threshold obtained as the eigen- 
vectors of 0(^x.y) with the largest growth rate and different 
angles of rotation of boundary against anisotropy direction a 
(value of abs{e) — \fl is plotted), (a)-(d) — anisotropy and 
boundaries are aligned (a = 0); (e)-(h) — the misalignment 
present (a = 7r/5); (i)-(l) — the misalignment has its largest 
possible value (a = 7r/4). (a),(e),(i) — x- component of the 
field polarization, (b),(f),(j) — its transverse Fourier trans- 
form; (c),(g),(k) — y- component of the field polarization, 
(d),(h),(l) — its transverse Fourier transform. The other rel- 
evant parameters are: I = 40/im 7^ = 0.01, 7p = 0, 5 = 30nm. 



B. Influence of inhomogeneities 

In this subsection, we consider an inffuence of pump 
inhomogeneities of different shape, amplitude and local- 
ization. As it was found earlier for devices with circular 
aperture [s^l, any slight inhomogeneities can change the 
sequence of a few first eigenfunctions in accordance with 
the order of their decay rates, preserving the shape of 
modes and their frequencies. In this connection, it is 
useful to investigate the subsequent eigenmodes of the 
linear eigenvalue problem 
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FIG. 2: The second mode of homogeneous cavity with pa- 
rameters as in Fig[T] (e)-(h). It can be made first by adding 
small inhomogeneity in current which is localized near 
right boundary of the device. 

The second eigenmode for the parameters of Fig.[lje)- 
(h) is shown in FigO The Fourier spectrum of this mode 
for both polarization components consists of four points: 
two pairs of two closed points. These pairs are situated 
approximately at the same places that two spots in Fig[T] 
(f), (h) and rotated by some angle against polarization di- 
rection. So, the main difference between two first modes 
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is the splitting of the Fourier harmonics leading to the 
strong modulation of stripes in the near field (Fig[2] (b), 

(d) ). The second eigenmode can be made leading one 
(one with the highest growth rate) by adding small inho- 
mogeneity to the current Sfi ~ 0.05(5(a;, y), where 6{x,y) 
is a probe function which is nonzero in the vicinity of the 
right border and zero everywhere else. 

If we decrease the detuning S between the gain maxima 
and the cavity resonance, the gain-loss dispersion mech- 
anism of Fourier harmonics selection becomes weaker 
[13, Ei and, as result, eigenmodes are more strongly de- 
termined by boundary conditions. The first mode for the 
large angle near tt/A does not have clear stripelike struc- 
tures in the center of device (Fig. [3] (a), (b)), and the sec- 
ond mode has the structure even less similar to stripes 
than previous one (Fig. [3] (c),(d)). However, stripelike 
structures do not disappear at all, as one can see in Fig. [3] 

(e) ,(f) for the next eigenmode. Moreover, the patterns 
with orthogonal direction of stripes can be easily exited 
(Fig. [21(g), (b)). This could be explained by the fact that 
these two families of stripes are degenerate when a = 0, 
and start to deviate from each other with increasing a. 
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FIG. 4: The subsequent eigenmodes for the same parameters 
as in Fig. (a)-(d) - 7th (e)-(h) - 21th, (i)-(l) - 26th eigen- 
mode. (a),(e),(i) — X- component of the field polarization, 
(b),(f),(j) — its transverse Fourier transform; (c),(g),(k) — 
y- component of the field polarization, (d),(h),(l) — its trans- 
verse Fourier transform. 



However, even in this case the transverse Fourier images 
of the field components look quite similar. This is con- 
nected to the fact, that structures in both polarizations 
belong to the same eigenmode and coupled to each other 
in the absence of strong inhomogeneities through bound- 
ary conditions (they both and some combination of their 
derivatives must vanish at the boundary), and through 
the Bragg refiectors. The exception is only when this 
connection is very weak (Fig. [T] (a)-(d)). 



FIG. 3: The first (a),(b); second (c,d); 5th and 6th eigen- 
modes of the device with the same parameters as in Fig. [1] 
but with smaller detuning (5 = lOnm), and angle of rota- 
tion of boundaries against anisotropy direction is a = 7r/4. 
(a),(c),(e),(g) - X component of the field, (b),(d),(f),(h) - y 
component of the field. 

Some subsequent eigenmodes with larger decay rate 
for the same parameters are presented in Fig. 3) It is 
evident that for small enough (but nonzero) anisotropy 
the modes can differ from each other in their symmetry 
properties essentially, as in the case of circular aperture 
[12| . One can find nearly unordered structures (Fig. 0] 
(e)-(h)) and the structures which can be considered as 
'scared' eigenmodes of quantum billiard (Fig. [5] (a)- (d), 
(i)-(l))- 

It is difficult to change the order of subsequent modes 
(especially to push one of them to the first position with 
highest grough rate) by introducing a small inhomogene- 
ity. Since stronger inhomogeneities needed for that mode 
structures can be changed also. 

It should be also noted that whereas patterns of the 
orthogonally-polarized component (weaker component) 
of the laser field are usually of the same shape, one can 
find some exceptions (compare Fig. (U^i) and Fig. (Hk)). 
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FIG. 5: The first two eigenmodes for the device with the same 
parameters as in Fig. [T] (a)-(d), but & = 8nm, and peak of 5/i 
in the aperture. (a),(e) — x- component of the field polariza- 
tion, (b),(f) — its transverse Fourier transform; (c),(g) — y- 
component of the field polarization, (d),(h) — its transverse 
Fourier transform. 

Up to now we have investigated the eigenmodes of ho- 
mogeneous devise, having in mind that the grough rate 
of a mode following the first one can become the largest 
one when a small perturbation of the current (5/i is added, 
which does not sufficiently change the shape of mode. 

Increasing the intensity of inhomogeneities leads to 
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changes in the shape of eigenmodes. As an example, 
we consider here strong inhomogeneities by introducing 
a peak or a hole into the laser aperture, that is quite 
natural for experiments. 

For the peak in current profile, the first cigenmode 
is defined by this spot only, as if the rest of aperture 
is empty (Fig. [5] (a)-(d)). But, the subsequent modes 
fill all the cavity as there is no peak in the aperture at 
all (Fig. [5] (e)-(h)). However, they are more disordered 
then in the absence of inhomogeneity. The corresponding 
eigenvalues are also more separated as it is described in 
the next section. 



(b) 




FIG. 6: The first eigenmode for the same parameters as in 
Fig. H] (a)-(d), but the boundary conditions (|C5|) are defined 
not on the boundaries of a square, but also on a circle inside 
the square (which creates the hole in the aperture). The hole 
is of 1/8 of the size of the laser and placed in the position 
(1/4, 1/4) in the whole aperture, (a) — x- component of the 
field polarization, (b) — its transverse Fourier transform; (c) 
— y- component of the field polarization, (d) — its transverse 
Fourier transform. 

On the other hand, if there is a strong hole in the 
current profile (which is for simplicity modeled by intro- 
ducing zero boundary conditions on the circle surround- 
ing the hole), the effect is not so visible, i.e. there is 
not a separation on the "hole" mode and the rest ones. 
Moveover, the shape of eigenmodes for x polarization is 
not changed so dramatically far from the hole (Fig. [H]), 
although the structure of weak component is strongly dis- 
ordered comparing to the nonperturbed case (Fig.[IJc)). 
However, the subsequent eigenmodes are more irregular 
that the first one even for the x-component. The problem 
becomes close to chaotic one that could be considered in 
the framework of wave chaos. 



IV. SIGNATURES OF QUANTUM CHAOS 

Generally speaking, the term 'quantum chaos' or 'wave 
chaos' is usually attributed to investigation of quantum 
systems, which posses chaotic features in classical limit 
|25j . In contrast to classical systems, they obey linear 
evolution equation (Schrodinger equation). So, the ob- 
ject of investigation is usually the set of eigenvalues and 
eigenf unctions (stationary states) of some linear opera- 
tor, describing the quantum system under consideration 
(usually Hamilton operator). At that, the very important 
role is played by eigenvalue distribution, which shows fin- 
gerprints of chaos in quantum systems when the corre- 
sponding classical systems are chaotic. 

Among the quantum systems, much attention is de- 
voted to the billiard systems [2^,[3^. The quantum coun- 



terpart of classical billiard is described by a wavefunction 
'(/'(x, y, t), which obeys Schrodinger equation: 



with boundary condition 



= o| 



dS ' 



(4) 



(5) 



on the boundary of the domain S. 

Stationary states of this problem are eigenfunctions of 
the transverse Laplace operator A = + 

iAiP = AV'. (6) 



Depending on the shape of domain 5, the correspond- 
ing classical problem can be either chaotic or regular. In 
the limit case of integrable system (for example, when the 
region has a rectangular shape), the corresponding eigen- 
value problem gives set of eigenvalues, obeying Poisson 
statistics of an eigenvalue spacing s: 



p{s) = exp{-s), 



(7) 



where p{s) is a probability for corresponding value of 
eigenfrequency distribution Si — ImXi — Im A^+i. 

For the opposite limit case of completely chaotic sys- 
tem, the corresponding distribution is the Wigner one: 



(8) 



For intermediate situations one can define other statis- 
tic famihes hmited by Eq. ^ and Eq. jS]) [H, HE HJ. 
One of fingerprints of quantum chaos in these intermedi- 
ate cases, as well as in the Eq. ([5]), that the maximum 
of p(s) is reached for s 7^ 0. This defines so called level 
repulsion phenomena for chaotic systems. 

The eigenfunctions of Hamilton operator also demon- 
strate fingerprints of chaos. In particular, one can ob- 
serve so called 'scarred' patterns [il, [13, [H, S SO], 
which are localized near (unstable) periodic trajecto- 
ries of corresponding classical systems. Among the 
scared patterns, there are other eigenfunctions, which are 
strongly irregular and fit all the area. 

It should be noted, that besides the quantum systems, 
the above mentioned framework is of the considerable 
interest for macro-systems, which are described at some 




or Eq. ([4]). Among 
W\ , microdisk lasers 



level of approximation by Ec 
them are microwave billiards j2£ 

[H,!!!] and vcsELs [H,!!!,!!! 

The eigenvalue statistics and other fingerprints of 
quantum chaos can give a criterion of a complexity of 
patterns near threshold. Indeed, if a system possess the 
Poisson statistics, the corresponding eigenfunctions are 
regular and stable against small perturbations, whereas 
for the Wigner one the eigenfunctions are much more 
irregular. 

In the present work, we describe the pattern forma- 
tion in VCSEL by linear operator O ([T|) obeying Eq. ((2) 
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This operator is considerably more complicated then the 
transverse Laplace operator from Eq. First of all, 
the eigenvalues of the operator O are complex. Their 
real parts describe the decay rate of the corresponding 
eigenf unctions, whereas the complex parts are their os- 
cillation frequencies determining an energy level distribu- 
tion. Examples of the eigenvalues distribution are pre- 
sented in Fig. [71 For homogeneous device, all eigenvalues 
are concentrated along one line and in average their de- 
cay rates are increased with the frequencies (Fig. [7] (a)). 
The closeness of decay rates of the adjacent modes allows 
to change their order by a small perturbation. It is worth 
noting that for the taken parameters the eigenvalues are 
grouped in clusters and become deviate with decreasing 
of the angle a. In the case of a local peak in the pump 
current the eigenvalues distribution confirms the separa- 
tion of the first and the rest modes (Fig. [7] (b)). 
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On the other hand, the operator ([T]) is anisotropic, with 
a preferential directions defined by anisotropy of VCSEL 
cavity. As one can see from (II|),(l2]), if the boundaries are 
defined parallel to x and y axes, one can separate vari- 
ables in the eigenvalue problem and therefore it becomes 
integrable. Hence, as one can see in Fig. [8l^a), when the 
boundaries and anisotropy direction coincides, the eigen- 
values obey more or less Poisson statistics. 

When these two directions (anisotropy and boundaries 
for square shape of the aperture) are not aligned, the 
variable separation is not possible anymore (as would 
be for the second order differential operator). As a re- 
sult, for a 7^ (Fig. [5{b)) one can clearly see that the 
maximum probability of the level spacing distribution 
is shifted to nonzero s, possessing a fingerprint of wave 
chaos. However, the system can not be considered as 
purely 'quantum chaotic' because the statistic is of inter- 
mediate type, neither purely Poisson nor Wigner one. It 
should be noted, that for isotropic system obeying Eq. ^ 
the eigenvalue spacing statistics remain Poisson for any 
values of a. In our case, the eigenvalue spacing has the 
maximum value again at zero for the angle a = 7r/4, that 
can be explained by possible degeneracy of the operator 
6 (Fig. He)). 

As a comparison the eigenlevel statistics is shown for 
boundary conditions of Fig. [6l It is clearly seen that 
for this case the statistics for spacing of imaginary parts 
of eigenvalues can be considered as Wigner one, which 
is quite expectable, because the same happens with 
the eigenvalues of the transverse Laplace operator. On 
the other hand, for the large peak in the aperture the 
eigenvalues do not obey the Wigner statistics anymore. 
For this case, the leading eigenvalues (eigenfunctions for 
which are shown in Fig.[5l[a)-(d) and can be described as 
eigenvalues of a disturbance only) are strongly separated 
from others (which are resembling Fig. [5Ke)-(h)), which 
obviously breaks the Wigner distribution. 

The other signatures of complexity of spatial structures 
can be found in fingerprints of quantum chaos demon- 
strated by the shape of eigenfunctions. As one can see, 
the neighboring eigenmodes for misaligned boundaries 
and anisotropy directions (a ^ 0) can differ sufficiently 
(Fig. [SI Fig. [4]) , including very irregular ones (Fig. 
(g)), which is a common property of systems possessing 
quantum chaos. On the other hand, beside complicated 
unordered structures, quantum chaos is characterized by 
'scared' eigenfunctions, which are located near unstable 
periodic trajectories of corresponding classical system. 
The example of eigenfunctions of operator O, resembling 
that kind of structures (Fig.[4|)(a)-(d),(i)-(l). Experimen- 
tally these structures were observed recently in VCSEL 



FIG. 7: (a) — Imaginary versus real part of few tens eigen- 
values \i of operator O having the largest value of Re \i for 
VCSEL with parameters as in Fig. [1] (e)-(h). (b) — Imagi- 
nary versus real part of about hundred eigenvalues \i of O 
for parameters close to ones used for Fig. [S] 



V. DISCUSSION AND CONCLUSION 

We have considered the pattern formation in wide 
aperture VCSEL near lasing threshold. Because of the 
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(a) 




FIG. 8: The statistics of approximately 100 first eigenmodes 
(with largest value of Re A) of VCSEL with parameters and 
boundaries as in Fig. [T] and Fig. [G] shown as stairstep graph, 
bold line. For a comparison Poisson statistics (Eq.O Wigner 
statistics (Eq. (Sjl is shown by dashed and dot-dashed lines, 
correspondingly, (a) corresponds to (a)-(d) in Fig.[T] i.e. a = 
0, (b) corresponds to (e)-(h) in Fig. [T] i.e. a = tt/5, (c) 
corresponds to (i)-(l) in Fig. [T] i.e. a — n/4, (d) corresponds 
to complicated region of Fig. [S] 



Bragg reflectors, the presence of polarization anisotropy 
creates a spatial anisotropy. This allows to investigate 
the structure of pattern near threshold by linearized or- 
der parameter equations, since even in a linear approxi- 
mation the mode with the smallest decay rate is nonde- 
generate (in contrast to spatially isotropic systems, where 
the whole family of modes has the same growth rate at 
threshold, and nonlinear competition is of strong impor- 
tance for the selection). Since our studies concentrated 
on a VCSEL with square aperture, the competition of 
boundary direction (which can be defined as a vector, 
parallel to one of the sides of the square), and spatial 
anisotropy is important mechanism for pattern selection. 

To measure the complexness of a pattern structure we 
use criteria from the topic of quantum chaos, where the 
main object of research is also an eigenvalue problem of 
some operator. One of fingerprints of complicated be- 
havior is an eigenvalue spacing statistics. 

In this work we have considered a statistics of imagi- 
nary parts of eigenvalues, describing frequencies of cor- 
responding eigenfunctions. We have shown that in the 
case of aligned directions of boundaries and spatial (or 
polarization) anisotropy, the above mentioned statistics 
is close to Poisson one. However, when the angle of rota- 
tion of spatial anisotropy against boundaries is not zero, 
the statistics is not Poisson anymore. In the later case, 
the statistics shows a signature of level spacing repulsion, 
i.e. maximum of eigenspacing distribution is not zero. 

The other fingenprints of quantum chaos can be found 
in the spatial shape of eigenfunctions, which can be very 
irregular, as well as 'scared-like' patterns. 

To compare these results with experiments, we also 
have investigated the influence of the inhomogeneities. 
The considered system is stable in the sense, that small 
perturbations only change the modes order preserving 
their shape, and not lead to dramatic change of the whole 
picture. In this sense, the consideration of subsequent 
modes (having higher decay rate) are also important, be- 
cause they can be considered as leading ones of the device 
perturbed by some (may be unknown) inhomogeneity. 

With increasing the inhomogeneity, the shape of eigen- 
modes also changes, as well as all the characteristics like 
the level spacing distribution. One can, for example, 
model the hole in the current by imposing zero boundary 
conditions somewhere on a circle located inside the squire 
aperture. In this case, the statistics evidently becomes 
Wigner, because the corresponding billiard is scattering 
one. On the other hand, large positive peak in the cur- 
rent leads to separation of eigenvalues and eigenfunctions 
into two classes, first of them belonging to 'peak family', 
with eigenfunctions localized near the peak, whereas the 
second family is distributed across the whole cavity aper- 
ture. 
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account propagation inside the complex VCSEL cavity. 
The material equations are based on spin- flip model [4l| , 
describing four level system with populations differences 
of relevant transitions determined by their sum N and 
difference n: 



APPENDIX A: BASIC SYSTEM 

The system under consideration was initially obtained 
in The equations for the field E = (i?^, Ey) take into 

I 



dE{xt,t) 



dt 



dN 

It 

dn 



-nAm + iVlE - TE - ikoE + k(1 + ia)GLn{A¥.), 
-N + fi- lm[{i - a)E*Ln{AE)], 
-j,d-Re[{i-a)WLniAE,)], 



(Al) 

(A2) 

(A3) 
(A4) 



r 



where 



A 



N in 

-in N 



and the operators M and G describes the modal losses 
and gain, correspondingly, VL describes the influence of 
the diffraction in the complex resonator of VCSEL. These 
operators take into account the influence of reflection 
from Bragg reflectors and are quite complicated. More 
details can be found in [ll|, [l^] ■ 
The operator 



Ln = l/\l 



d-n 

7 



describes the Lorentz shape of gain contour with a detun- 
ing S between the gain maximum and cavity resonance 
frequency, k is the field decay rate in the VCSEL's cavity. 

The operator T describes the inner anisotropy of VC- 
SEL's cavity: 

exp(7Q + ijp) 

exp(-(7a + ijp)) 

where ja and 7^ is an amplitude and phase anisotropy, 
respectively. 

APPENDIX B: VECTORIAL EIGENVALUE 
PROBLEM DERIVATION 

To derive the linear evolution equation in the form, 
dE 
'dt 



(Bl) 



(which is then reduced to Eq.[5]) using the basic equations 
(lAlllASI) . we take into account that at the laser thresh- 
old two branches of the steady state solutions cross each 
other: the zero solution (E = 0), and the nonzero lasing 
one. Because of the lasing solution at the cross-section 
point is characterized by E = and N = the further 
analysis is drastically simplified giving the following di- 
agonal operator of the linearized problem: 



T 



On 
622 



(B2) 



where each On is in turn the matrix operator: 



Aiii) 
"-^11 

o. 



) 



(B3) 



The diagonal form of (jB2j) shows, that the field E and 
the variables N ,n are independent of each other at the 
laser threshold, and the spectrum for O22 lies entirely in 
the half-plane Re A < 0, hence we can consider only the 
equation for the field Eq. IB II with O = On. O includes 
the current density at threshold N — ^. From (jAip . (|B2p . 
the explicit value for the operator O is defined by the 
following action on the field E: 
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O — ^kM + if7 — r — iKa + k(1 + ia)fiGL(i. 

I 



However, the expression (jB4|) is still nonlocal, because 
all the operators in (jB4[) are integrodifferential. The 
equation can be made local in the manner of [12]. The 
procedure used for a scalar operator in [12], must be ap- 
plied to each of the operators Oj-j^"*, which are the scalar 

components of O. Because the diagonal part of (jB3[) 
plays the main role, it is approximated up to the fourth 
order of kj_, whereas the operators O^-j^'* for i ^ j are 
approximated with terms of the second order of kj^ . 

In the presence of index inhomogeneities, Ct = 
^homogen + HSn (whcrc I = TLo/uo IS cxprcsscd via the 
optical frequency lu, round trip time of the cavity rand 
the mean index uq. The pump current inhomogeneities 
have to be included in the coefficient of the last term 

M — f^homogen ~^ S^. 



APPENDIX C: NUMERICAL 
IMPLEMENTATION OF EIGENVALUE 
PROBLEM 

To solve the problem, the Matlab PDE toolbox imple- 
mentation ^2] of Arnoldi method for the system of par- 
tial differential equations of the second order has been 
used. For that, the operators of the fourth order O^l^^ 
have been represented via multiplication of two operators 
Pii, Pi2 of the second order with constant coefficients: 

^'ti^^ ~ PilPi2 ~ hi — Pi2Pil — hi (CI) 

where Pn^W ■ {c{l «) V), P,2 = V ■ {c^i ® 
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